Code to load tools and data

# TOOLS 
  require(here)
  source(here::here('R/tools.R'))
  
  require(arm) 
  require(effects)
  require(ggpubr)
  require(ggsci) 
  require(gridExtra)
  require(magrittr)
  require(multcomp)
  require(rptR) 
  require(stringi)
  require(viridis)


  # constants
    round_ = 3 # number of decimal places to round model coefficients
    nsim = 5000 # number of simulations to extract estimates and 95%CrI
    ax_lines = "grey60" # defines color of the axis lines
    colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
  
  # functions
    getime = function (x) {ifelse(is.na(x), as.numeric(NA), as.numeric(difftime(x, trunc(x,"day"), units = "hours")))}
    
    getDay = function (x) {as.Date(trunc(x, "day"))}
   
    # for R output based on sim
      R_out = function(name = "define", model = m, nsim = 5000){
       bsim <- sim(model, n.sim=nsim)  
       l=data.frame(summary(model)$varcor)
       l = l[is.na(l$var2),]
       l$var1 = ifelse(is.na(l$var1),"",l$var1)
       l$pred = paste(l$grp,l$var1)

       q50={}
       q025={}
       q975={}
       pred={}
       
       # variance of random effects
       for (ran in names(bsim@ranef)) {
         #ran =names(bsim@ranef)[1]
         ran_type = l$var1[l$grp == ran]
         for(i in ran_type){
            # i = ran_type[2]
          q50=c(q50,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.5)))
          q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.025)))
          q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.975)))
          pred= c(pred,paste(ran, i))
          }
         }
       # residual variance
       q50=c(q50,quantile(bsim@sigma^2, prob=c(0.5)))
       q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
       q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
       pred= c(pred,'Residual')

       ci = c(round(100*q025/sum(q025))[1], round(100*q975/sum(q975))[1])
       ci = ci[order(ci)]
       
       ri=data.table(model = name, repeatability=paste0(round(100*q50/sum(q50)),'%')[1], CI = paste0(paste(ci[1], ci[2], sep ="-"), '%'))
       
       
       return(ri)
       }

# DATA  
  d = data.table(read_excel(here::here('Data/motility.xlsx'), sheet = 1))
  s = data.table(read_excel(here::here('Data/sampling_2021_cleaned.xlsx'), sheet = 1))
  s = s[!is.na(recording)]

  # add morph and age
    m = data.table(read_excel(here::here('Data/ruff_males_Seewiesen.xlsx'), sheet = 1))#, range = "A1:G161"))
    m[, hatch_year:=as.numeric(substr(Hatch_date,1,4)) ]
    m[, age := 2021-hatch_year]

    m = m[Ind_ID %in%unique(s$DB_ID[!is.na(s$DB_ID)]),.(Ind_ID, Morph, age)]
    s = merge(s,m, by.x = 'DB_ID', by.y = 'Ind_ID', all.x = TRUE)
    d = merge(d, s[,.(bird_ID, Morph, age)], by.x = 'ID', by.y = 'bird_ID', all.x = TRUE)
    d = (d[!duplicated(d)]) # shouldn't be necessary, but for some reason the above call duplicates the dataset
    
    d[is.na(Morph), Morph := 'Zebra finch']
    d[Morph == 'F', Morph := 'Faeder']
    d[Morph == 'I', Morph := 'Independent']
    d[Morph == 'S', Morph := 'Satellite']
    d[is.na(issues), issues := 'zero']

    #nrow(d) # N 139

  # prepare for correlations and repeatability
    dr = d[ID%in%ID[duplicated(ID)]]
    drw = reshape(dr[,.(date,ID,VAP,VSL,VCL, motileCount, Morph, age)], idvar = c('ID','Morph','age'), timevar = 'date', direction = "wide")  

Sample sizes

    #table(d$ID, d$date)

    # June
    j = s[datetime>as.POSIXct('2021-05-15')]
    nrow(j) # 108 recordings
    length(unique(j$bird_ID)) # from 99 sampled individuals

    nrow(d[date == 'June']) # 93 recordings analyzed
    length(unique(d[date == 'June', ID])) # for 93 individuals

    jj = unique(j$bird_ID)
    dd = unique(d[date == 'June', ID])
    jj[!jj%in%dd]
    dd[!dd%in%jj] # all IDs in motility measurements are from the sampled 

    # May
    m = s[!datetime>as.POSIXct('2021-05-15')]
    nrow(m) # 53 recordings
    length(unique(m$bird_ID)) # 47 of 
    nrow(d[date == 'May']) # 46 + 1 (A01718) without tracking data

Exploration

Densities

Velocity

    g1 = 
    ggplot(d,aes(x = VCL, col = Morph)) + 
      geom_density() + 
      xlab('Velocity μm/s') + 
      xlim(c(0,72)) + 
      ggtitle('Curvilinear') +
      theme_bw() +
      theme(
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        legend.position = c(0.15, 0.65),
        legend.text=element_text(size=6),
        legend.title=element_text(size=7),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
    g2 = 
    ggplot(d,aes(x = VAP, col = Morph)) + 
      geom_density() + 
      xlim(c(0,72)) +
      ggtitle('Average path') + 
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_text(hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.text.x = element_blank())
    
    g3 =  
    ggplot(d,aes(x =  VSL, col = Morph)) + 
      geom_density() + 
      xlim(c(0,72)) +
      xlab('Velocity μm/s') + 
      ggtitle('Straight line') +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_text(color = 'white', hjust = 0.5)
          )

    grid.draw(rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-density.png'),rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"), width = 7*1.5, height =8*1.5, units = 'cm')    

Velocity - date

    g1 = 
    ggplot(d,aes(x = VCL, col = Morph)) + 
      geom_density() + 
      xlab('Velocity μm/s') + 
      xlim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~date, ncol = 2) +
      theme_bw() +
      theme(
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        legend.position = c(0.15, 0.55),
        legend.text=element_text(size=6),
        legend.title=element_blank(),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
    g2 = 
    ggplot(d,aes(x = VAP, col = Morph)) + 
      geom_density() + 
      xlim(c(0,72)) +
      ggtitle('Average-path') + 
       facet_wrap(~date, ncol = 2) +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_text(hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.text.x = element_blank())
    
    g3 =  
    ggplot(d,aes(x =  VSL, col = Morph)) + 
      geom_density() + 
      xlim(c(0,72)) +
      xlab('Velocity μm/s') + 
      ggtitle('Straight-line') +
      facet_wrap(~date, ncol = 2) +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_text(color = 'white', hjust = 0.5)
          )

    grid.draw(rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-density-date.png'),rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"), width = 7*1.5, height =8*1.5, units = 'cm')  

N tracked sperm

    g1 =      
    ggplot(d,aes(x = motileCount, col = Morph)) + 
      geom_density() + 
      xlab('N tracked sperm cells') + 
      theme_bw() +
      theme(
        axis.title.y = element_blank(),
        legend.position = c(0.85, 0.75),
        legend.text=element_text(size=6),
        legend.title=element_text(size=7),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
    
    #densityplot(~motileCount, group = date, data = d, auto.key = TRUE)
    
    g2 = 
    ggplot(d,aes(x = motileCount, col = date)) + 
      geom_density() + 
      xlab('N tracked sperm cells') + 
      scale_colour_discrete(guide = guide_legend(reverse = TRUE)) + 
      theme_bw() +
      theme(
        axis.title.y = element_blank(),
        legend.position = c(0.85, 0.75),
        legend.text=element_text(size=6),
        legend.title=element_text(size=7),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
    
    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), size = "first"))

    ggsave(here::here('Output/Motility-density_N.png'), cbind(ggplotGrob(g1), ggplotGrob(g2), size = "first"), width = 8*1.5, height =4*1.5, units = 'cm')    
    
    summary(d$motileCount) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     5.0    87.5   182.0   202.2   299.5   562.0

Age

    g1 = 
    ggplot(d,aes(x = age, col = Morph)) + 
      geom_density() + 
      xlab('Age') + 
      theme_bw() +
      theme(
        legend.position = c(0.85, 0.85),
        legend.text=element_text(size=6),
        legend.title=element_text(size=7),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
     
   g2 = 
    ggplot(d,aes(x = age, col = date)) + 
      geom_density() + 
      xlab('Age') + 
      scale_colour_discrete(guide = guide_legend(reverse = TRUE)) + 
      theme_bw() +
      theme(
        axis.title.y = element_blank(),
        legend.position = c(0.85, 0.85),
        legend.text=element_text(size=6),
        legend.title=element_text(size=7),
        legend.key = element_rect(colour = NA, fill = NA),
        legend.key.height= unit(0.5,"line"),
        legend.key.width = unit(0.25, "cm"),
        legend.margin = margin(0,0,0,0, unit="cm"),
        legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
        legend.background = element_blank(),
        plot.title = element_text(size=9, hjust = 0.5))
      
    g3 = 
     ggplot(d,aes(x = age, col = Morph)) + 
      geom_density() + 
      xlab('Age (log)') + 
      theme_bw() +
      scale_x_continuous(trans = 'log') + 
      theme(
        axis.title.y = element_blank(),
        legend.position = "none",
        plot.title = element_text(size=9, hjust = 0.5))
    
    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g3), ggplotGrob(g2), size = "first"))

    ggsave(here::here('Output/Motility-density_age.png'), cbind(ggplotGrob(g1), ggplotGrob(g3),ggplotGrob(g2), size = "first"), width = 10*1.5, height =4*1.5, units = 'cm') 

Histograms

Velocity

    g1 = 
    ggplot(d,aes(x = VCL, col = Morph)) + 
      geom_histogram() + 
      coord_cartesian(xlim = c(0,72), ylim = c(0,20)) +
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
         axis.title.x = element_text(colour= 'white', hjust = 0.75),
         plot.title = element_text(size=9, hjust = 0.5))
    g2 = 
    ggplot(d,aes(x = VAP, col = Morph)) + 
      geom_histogram() + 
       coord_cartesian(xlim = c(0,72), ylim = c(0,20)) +
      xlab('Velocity μm/s') + 
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.text.y = element_blank(),
          axis.title.y = element_blank()
          )
    
    g3 =  
    ggplot(d,aes(x =  VSL, col = Morph)) + 
      geom_histogram() + 
      coord_cartesian(xlim = c(0,72), ylim = c(0,20)) +
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank()
          )

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-hist.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')       

N and age

    g1 = ggplot(d,aes(x = motileCount, col = Morph)) + 
      geom_histogram() + 
      #coord_cartesian(xlim = c(0,72), ylim = c(0,20)) +
      xlab('N tracked sperm cells') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none"
         )

    g2 = ggplot(d,aes(x = age, col = Morph)) + 
      geom_histogram() + 
      #coord_cartesian(xlim = c(0,72), ylim = c(0,20)) +
      xlab('Ages') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none"
         )  
    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), size = "first"))

    ggsave(here::here('Output/Motility-hist_N&Age.png'), cbind(ggplotGrob(g1), ggplotGrob(g2), size = "first"), width = 7*1.5, height =7*1.5, units = 'cm')  

Correlations

Velocity ~ N all

    g1 = 
    ggplot(d,aes(x = motileCount, y = VCL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),formula = y ~ poly(x,2), size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
         axis.title.x = element_text(color="white"), 
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = motileCount, y = VAP)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),formula = y ~ poly(x,2), size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = motileCount, y = VSL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),formula = y ~ poly(x,2), size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corSpermN.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  

Velocity ~ N by sampling date

    g1 = 
    ggplot(d,aes(x = motileCount, y = VCL, col = date)) + 
      geom_point(pch = 21) + 
      stat_smooth(method = "lm", formula = y ~ poly(x,2), size = 1, se = FALSE)+
      #stat_cor(method="pearson",size = 2) +
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.x = element_text(color = "white"),
        plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = motileCount, y = VAP, col = date)) + 
      geom_point(pch = 21) + 
      stat_smooth(method = "lm", formula = y ~ poly(x,2), size = 1, se = FALSE)+
      #stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = motileCount, y = VSL, col = date)) + 
      geom_point(pch = 21) + 
      stat_smooth(method = "lm", formula = y ~ poly(x,2), size = 1, se = FALSE) +
      #stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corSpermN_by_Date.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  

Velocity ~ log(N)

    g1 = 
    ggplot(d,aes(x = motileCount, y = VCL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      scale_x_continuous(trans='log10')+
      ylab('Velocity μm/s') + 
      ylim(c(0,70)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = motileCount, y = VAP)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      scale_x_continuous(trans='log10')+
      ylim(c(0,70)) +
      ggtitle('Average path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = motileCount, y = VSL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph),size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      scale_x_continuous(trans='log10')+
      ylim(c(0,70)) + 
      ggtitle('Straight line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corlogSpermN.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')    

Velocity ~ date

    g1 = 
    ggplot(d,aes(x = date, y = VCL, col = Morph)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.x = element_text(color = "white"), 
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = date, y = VAP)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = date, y = VSL, col = Morph)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corDate.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  

Velocity ~ age

    g1 = 
    ggplot(d[date == 'June'],aes(x = age, y = VCL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), formula = y ~ poly(x,2), size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
         axis.title.x = element_text(color = "white"), 
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d[date == 'June'],aes(x = age, y = VAP)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), formula = y ~ poly(x,2), size = 1, se = FALSE)+
      stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d[date == 'June'],aes(x = age, y = VSL)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), formula = y ~ poly(x,2), size = 1, se = FALSE) +
      stat_cor(method="pearson",size = 2) +
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corage.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  

Velocity ~ issues - Morph

    g1 = 
    ggplot(d,aes(x = issues, y = VCL, col = Morph)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      facet_wrap(~Morph, ncol = 1)  +
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.x = element_text(color = "white"), 
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = issues, y = VAP)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      facet_wrap(~Morph, ncol = 1)  +
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = issues, y = VSL, col = Morph)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 2)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      facet_wrap(~Morph, ncol = 1)  +
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corIssues_Morph.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  

Velocity ~ issues

    g1 = 
    ggplot(d,aes(x = issues, y = VCL, col = issues)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), fill ='grey40', dotsize = 1)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      ylab('Velocity μm/s') + 
      ylim(c(0,72)) + 
      ggtitle('Curvilinear') +
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.x = element_text(color = "white"), 
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d,aes(x = issues, y = VAP, col = issues)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), fill ='grey40', dotsize = 1)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      ylim(c(0,72)) +
      ggtitle('Average-path') + 
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      xlab('Issues during sperm tracking')+
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d,aes(x = issues, y = VSL, col = issues)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(),  fill ='grey40', dotsize = 1)+
      #geom_boxplot(col = 'grey',alpha = 0.2,position = position_dodge())+
      #geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2,position = position_dodge())+
      ylim(c(0,72)) + 
      ggtitle('Straight-line') +
      guides(x =  guide_axis(angle = -45)) +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-corIssues.png'),cbind(gg1,gg2,gg3, size = "first"), width = 18, height =5, units = 'cm')  

N ~ date

    ggplot(d,aes(x = date, y = motileCount, col = Morph)) +  
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph), fill ='grey40', dotsize = 1)+
      geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) + 
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      scale_colour_viridis(discrete=TRUE)+
      theme_bw() 

    ggsave(here::here('Output/Motility-corSpermN_Date.png'), width = 13, height =6, units = 'cm')  

N ~ age

    g1 = 
    ggplot(d[date == 'May' & Morph!='Zebra finch'], aes(x = age, y = motileCount)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), size = 1) +
      stat_cor(method="pearson",size = 2) +
      ylab('N sperm cell tracked') + 
      coord_cartesian(xlim = c(1,13), ylim = c(0,600)) + 
      ggtitle('May') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
         axis.title.x = element_text(color = "white"),
         plot.title = element_text(size=9, hjust = 0.5))
    
    g2 = 
    ggplot(d[date == 'June' & Morph!='Zebra finch'],aes(x = age, y = motileCount)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), size = 1) +
      stat_cor(method="pearson",size = 2) +
      coord_cartesian(xlim = c(1,13), ylim = c(0,600)) + 
      ggtitle('June') + 
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank())
    
    g3 =  
    ggplot(d[Morph!='Zebra finch'],aes(x = age, y = motileCount)) + 
      geom_point(pch = 21, aes(col = Morph)) + 
      stat_smooth(method = "lm", aes(col = Morph), size = 1) +
      stat_cor(method="pearson",size = 2) +
      coord_cartesian(xlim = c(1,13), ylim = c(0,600)) + 
      ggtitle('Both') +
      facet_wrap(~Morph, ncol = 1)  +
      theme_bw() +
      theme(legend.position = "none",
          plot.title = element_text(size=9, hjust = 0.5),
          axis.title.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank()
          )

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility_corN-age.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =8*1.5, units = 'cm')  


Correlations between May & June

    g1 = 
    ggplot(drw, aes(x = VCL.May, y = VCL.June)) +
    facet_wrap(~Morph, ncol = 1)  +
    stat_smooth(method = 'lm', aes(col = Morph))+geom_point(pch = 21, aes(col = Morph))+
    stat_cor(method="pearson",size = 2) +
    geom_abline(b = 1, col = 'red', lty = 3) + 
    xlim(c(min(c(drw$VCL.May,drw$VCL.June)), max(c(drw$VCL.May,drw$VCL.June)))) +  ylim(c(min(c(drw$VCL.May,drw$VCL.June)), max(c(drw$VCL.May,drw$VCL.June)))) + 
    ggtitle('Curvilinear')+
    ylab('Velocity in June [μm/s[') + 
    theme_bw() + 
    theme(legend.position = "none",
        axis.text = element_text(size=7), 
        axis.title.x = element_text(size = 8, color = 'white'),
        axis.title.y = element_text(size = 8),
        strip.text.x = element_text(size = 7),
        plot.title = element_text(size=8, hjust = 0.5))
    
    g2 = 
    ggplot(drw, aes(x = VAP.May, y = VAP.June)) +
    facet_wrap(~Morph, ncol = 1)  +
    stat_smooth(method = 'lm', aes(col = Morph))+geom_point(pch = 21, aes(col = Morph))+
    stat_cor(method="pearson",size = 2) +
    geom_abline(b = 1, col = 'red', lty = 3) + 
    xlim(c(min(c(drw$VAP.May,drw$VAP.June)), max(c(drw$VAP.May,drw$VAP.June)))) +  ylim(c(min(c(drw$VAP.May,drw$VAP.June)), max(c(drw$VAP.May,drw$VAP.June)))) + 
    ggtitle('Average path')+
    xlab('Velocity in May [μm/s]') + 
    theme_bw() + 
    theme(legend.position = "none",
        axis.text = element_text(size=7), 
        axis.title.y = element_blank(), 
        axis.title.x = element_text(size = 8, hjust = 0.5),
        strip.text.x = element_text(size = 7),
        plot.title = element_text(size=8, hjust = 0.5))

    g3 = 
    ggplot(drw, aes(x = VSL.May, y = VSL.June)) +
    facet_wrap(~Morph, ncol = 1)  +
    stat_smooth(method = 'lm', aes(col = Morph))+geom_point(pch = 21, aes(col = Morph))+
    stat_cor(method="pearson",size = 2) +
    geom_abline(b = 1, col = 'red', lty = 3) + 
    xlim(c(min(c(drw$VSL.May,drw$VSL.June)), max(c(drw$VSL.May,drw$VSL.June)))) +  ylim(c(min(c(drw$VSL.May,drw$VSL.June)), max(c(drw$VSL.May,drw$VSL.June)))) + 
    ggtitle('Straight line')+
    theme_bw() + 
    theme(legend.position = "none",
      axis.text = element_text(size=7), 
      axis.title.x = element_blank(), 
      axis.title.y = element_blank(), 
      strip.text.x = element_text(size = 7),
      plot.title = element_text(size=8, hjust = 0.5))

    grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "first"))

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility-cor_MayJune.png'),cbind(gg1,gg2,gg3, size = "first"), width = 7*1.5, height =7*1.5, units = 'cm')  

Repeatability

    # prepare
      velo = 'VAP'
        m = lmer(VAP ~ 1+(1|ID), dr)
        Rf_VAP = R_out(velo)
        Rf_VAP[, method_CI:='arm package']
        names(Rf_VAP)[1] = 'velocity'
        Rf_VAP[, pred:= as.numeric(substr(repeatability,1,nchar(repeatability)-1))]
        Rf_VAP[, lwr:= as.numeric(substr(CI,1,2))]
        Rf_VAP[, upr:= as.numeric(substr(CI,4,5))]

        R = rpt(VAP ~ (1 | ID), grname = "ID", data = dr, datatype = "Gaussian")#, nboot = 0, npermut = 0)
        RR = data.table(merge(data.frame(velocity =velo), paste0(round(R$R*100),'%'))) %>% setnames(new = c('velocity', 'repeatability'))
        RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')] 
        RR[, method_CI := 'rpt package']
        RR[, pred := 100*R$R]
        RR[, lwr := 100*R$CI_emp[1]]
        RR[, upr := 100*R$CI_emp[2]]
        RR_VAP = RR

        VAP = rbind(Rf_VAP, RR_VAP)

      velo = 'VSL'
        m = lmer(VSL ~ 1+(1|ID), dr)
        Rf_VSL = R_out(velo)
        Rf_VSL[, method_CI:='arm package']
        names(Rf_VSL)[1] = 'velocity'
        Rf_VSL[, pred:= as.numeric(substr(repeatability,1,nchar(repeatability)-1))]
        Rf_VSL[, lwr:= as.numeric(substr(CI,1,2))]
        Rf_VSL[, upr:= as.numeric(substr(CI,4,5))]

        R = rpt(VSL ~ (1 | ID), grname = "ID", data = dr, datatype = "Gaussian")#, nboot = 0, npermut = 0)
        RR = data.table(merge(data.frame(velocity =velo), paste0(round(R$R*100),'%'))) %>% setnames(new = c('velocity', 'repeatability'))
        RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')] 
        RR[, method_CI := 'rpt package']
        RR[, pred := 100*R$R]
        RR[, lwr := 100*R$CI_emp[1]]
        RR[, upr := 100*R$CI_emp[2]]
        RR_VSL = RR

        VSL = rbind(Rf_VSL, RR_VSL)

      velo = 'VCL'
        m = lmer(VCL ~ 1+(1|ID), dr)
        Rf_VCL = R_out(velo)
        Rf_VCL[, method_CI:='arm package']
        names(Rf_VCL)[1] = 'velocity'
        Rf_VCL[, pred:= as.numeric(substr(repeatability,1,nchar(repeatability)-1))]
        Rf_VCL[, lwr:= as.numeric(substr(CI,1,2))]
        Rf_VCL[, upr:= as.numeric(substr(CI,4,5))]

        R = rpt(VCL ~ (1 | ID), grname = "ID", data = dr, datatype = "Gaussian")#, nboot = 0, npermut = 0)
        RR = data.table(merge(data.frame(velocity =velo), paste0(round(R$R*100),'%'))) %>% setnames(new = c('velocity', 'repeatability'))
        RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')] 
        RR[, method_CI := 'rpt package']
        RR[, pred := 100*R$R]
        RR[, lwr := 100*R$CI_emp[1]]
        RR[, upr := 100*R$CI_emp[2]]
        RR_VCL = RR

        VCL = rbind(Rf_VCL, RR_VCL)

      r = rbind(VCL,VAP,VSL) 

      r[velocity == 'VCL', velocity := 'Curvilinear']
      r[velocity == 'VAP', velocity := 'Average path']
      r[velocity == 'VSL', velocity := 'Straight line']
    # plot    
      ggplot(r, aes(x = velocity, y = pred, col = method_CI)) +
            geom_errorbar(aes(ymin = lwr, ymax = upr, col = method_CI), width = 0, position = position_dodge(width = 0.4) ) +
            #ggtitle ("Sim based")+
            geom_point(position = position_dodge(width = 0.4)) +
            #scale_fill_brewer(palette = "Set1", guide = guide_legend(reverse = TRUE)) +
            #scale_color_brewer(palette = "Set1", guide = guide_legend(reverse = TRUE))  +
            scale_color_viridis(discrete=TRUE, guide = guide_legend(reverse = TRUE))  +
            scale_fill_viridis(discrete=TRUE, guide = guide_legend(reverse = TRUE)) + 
            #scale_shape(guide = guide_legend(reverse = TRUE)) + 
            #geom_hline(yintercept = 99, col = 'red')+
            labs(x = 'Velocity', y = "Repeatability [%]")+
            scale_y_continuous(expand = c(0, 0), lim = c(0,80)) +
            #scale_y_continuous(breaks = c(65, 70, 80, 85, 90, 95, 100))+
            #geom_text(y =99, x =0.6, label = '99', col = 'red', size = 3)+
            coord_flip()+
            theme_bw() +
            theme(
                #axis.title.x = element_blank(),
                legend.text=element_text(size=6),
                legend.title=element_text(size=7),
                legend.key = element_rect(colour = NA, fill = NA),
                legend.key.height= unit(0.5,"line"),
                legend.key.width = unit(0.25, "cm"),
                legend.margin = margin(0,0,0,0, unit="cm"),
                legend.box.margin = margin(l = -6), #legend.justification = c(-1,0),
                legend.background = element_blank()
                )

      ggsave(here::here('Output/Motility-repeatability.png'), width = 7*1.5, height =2.75*1.5, units = 'cm')

    # table  
      knitr::kable(r[,1:4])
velocity repeatability CI method_CI
Curvilinear 47% 45-48% arm package
Curvilinear 48% 19-69% rpt package
Average path 35% 32-37% arm package
Average path 36% 7-60% rpt package
Straight line 24% 21-26% arm package
Straight line 24% 0-51% rpt package


Analyses of interest

Compare May, June, May & June

  # plot May   
    g1= 
    ggplot(d[date == 'May'],aes(x = Morph, y = VCL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(20,72)) +
        ylab('VCL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        ggtitle("May")+
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )

    g2= 
    ggplot(d[date == 'May'],aes(x = Morph, y = VSL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(8,47)) +
        ylab('VSL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )
    
    g3= 
    ggplot(d[date == 'May'],aes(x = Morph, y = VAP)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(13,51)) +
        ylab('VAP [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_text(color = "white"), 
          #axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )   

    grid.draw(rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"))
        #grid.arrange(g1,g2)
    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    #ggsave(here::here('Output/Motility_Morp_boxplots_June.png'),rbind(gg1,gg2,gg3, size = "last"), width = 7*1.5, height =15*1.5, units = 'cm')  
    #summary(glht(m)) 
    #plot(allEffects(m))
  # plot June   
    dj = d[date == 'June']
    # dummy row to include zebra finch in the plot
    df = d[date == 'May' & Morph == 'Zebra finch'][1,]
    df$date = 'June'
    df$VAP = 100
    df$VCL = 100
    df$VSL = 100

    dj = rbind(dj,df)

    g4= 
    ggplot(dj,aes(x = Morph, y = VCL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 0.85)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(20,72)) +
        ylab('VCL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        ggtitle("June")+
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),  
          plot.title = element_text(size=9)
          )

    g5= 
    ggplot(dj,aes(x = Morph, y = VSL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 0.45)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(8,47)) +
        ylab('VSL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),  
          plot.title = element_text(size=9)
          )
    
    g6= 
    ggplot(dj,aes(x = Morph, y = VAP)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 0.45)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(13,51)) +
        ylab('VAP [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          #axis.title.x = element_blank(), 
          #axis.text.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),  
          plot.title = element_text(size=9)
          )   

    
    grid.draw(rbind(ggplotGrob(g4), ggplotGrob(g5), ggplotGrob(g6), size = "last"))
        #grid.arrange(g1,g2)
    gg4 <- ggplotGrob(g4)
    gg5 <- ggplotGrob(g5) 
    gg6 <- ggplotGrob(g6) 
    #ggsave(here::here('Output/Motility_Morp_boxplots_May.png'),rbind(gg4,gg5,gg6, size = "last"), width = 7*1.5, height =15*1.5, units = 'cm')  
  # plot ALL  together 
    g7= 
    ggplot(d,aes(x = Morph, y = VCL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(20,72)) +
        ylab('VCL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        ggtitle("May & June")+
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(), 
          plot.title = element_text(size=9)
          )

    g8= 
    ggplot(d,aes(x = Morph, y = VSL)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(8,47)) +
        ylab('VSL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          axis.title.y = element_blank(), 
          axis.text.y = element_blank(), 
          plot.title = element_text(size=9)
          )
    
    g9= 
    ggplot(d,aes(x = Morph, y = VAP)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
        geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
        stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_viridis(discrete=TRUE)+
        coord_cartesian(ylim = c(13,51)) +
        ylab('VAP [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_text(color = "white"), 
          #axis.title.y = element_text(color = "white"), 
          axis.title.y = element_blank(), 
          axis.text.y = element_blank(), 
          #axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )   

    grid.draw(rbind(ggplotGrob(g7), ggplotGrob(g8), ggplotGrob(g9), size = "last"))
        #grid.arrange(g1,g2)
    gg7 <- ggplotGrob(g7)
    gg8 <- ggplotGrob(g8) 
    gg9 <- ggplotGrob(g9) 
    #ggsave(here::here('Output/Motility_Morp_boxplots.png'),rbind(gg1,gg2,gg3, size = "last"), width = 7*1.5, height =15*1.5, units = 'cm')
  # combine
    may = rbind(gg1,gg2,gg3, size = "last")
    june = rbind(gg4,gg5,gg6, size = "last")
    all = rbind(gg7,gg8,gg9, size = "last")

    grid.draw(cbind(may,june,all, size = "first"))

    ggsave(here::here('Output/Motility_Morp_boxplots_MayJuneAll.png'),cbind(may,june,all, size = "last"), width = 10*1.5, height =10*1.5, units = 'cm')

Compare May & June

    g1= 
    ggplot(d,aes(x = Morph, y = VCL, col = date)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     dotsize = 1, aes(col = date, fill = date),
                     position = position_dodge(width = 0.3))+
        geom_boxplot(outlier.color = NA, fill = NA, alpha = 0.2, width = 0.3, position = position_dodge(1.1)) + 
        #stat_summary(fun=mean, geom="point", color="red", fill="red") +
        #scale_color_viridis(discrete=TRUE)+
        scale_fill_npg( alpha = 0.2)+
        scale_color_npg()+
        #scale_color_viridis(discrete=TRUE)+
        #scale_fill_viridis(discrete=TRUE, alpha = 0.2)+
        ylab('VCL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "top",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )

    g2= 
    ggplot(d,aes(x = Morph, y = VSL, col = date)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     dotsize = 1, aes(col = date,fill = date),
                     position = position_dodge(width = 0.3))+
        geom_boxplot(outlier.color = NA, fill = NA, alpha = 0.2, width = 0.3, position = position_dodge(1.1)) + 
        #stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_npg( alpha = 0.2)+
        scale_color_npg()+
        #scale_color_viridis(discrete=TRUE)+
        #scale_fill_viridis(discrete=TRUE, alpha = 0.2)+
        ylab('VSL [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          axis.title.x = element_blank(), 
          axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )
    
    g3= 
    ggplot(d,aes(x = Morph, y = VAP, col = date)) + 
        geom_dotplot(binaxis = 'y', stackdir = 'center',
                     dotsize = 1, aes(col = date,fill = date),
                     position = position_dodge(width = 0.3))+
        geom_boxplot(outlier.color = NA, fill = NA, alpha = 0.2, width = 0.3, position = position_dodge(1.1)) + 
        #stat_summary(fun=mean, geom="point", color="red", fill="red") +
        scale_fill_npg( alpha = 0.2)+
        scale_color_npg()+
        #scale_color_viridis(discrete=TRUE)+
        #scale_fill_viridis(discrete=TRUE, alpha = 0.2)+
        ylab('VAP [μm/s]') +
        theme_bw() +
        guides(x =  guide_axis(angle = -45)) +
        theme(
          legend.position = "none",
          #axis.title.x = element_blank(), 
          #axis.text.x = element_blank(),
          plot.title = element_text(size=9)
          )   

    grid.draw(rbind(ggplotGrob(g1), ggplotGrob(g2), ggplotGrob(g3), size = "last"))

      #grid.arrange(g1,g2)

    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    ggsave(here::here('Output/Motility_Morp_boxplots_date.png'),rbind(gg1,gg2,gg3, size = "last"), width = 7*1.5, height =15*1.5, units = 'cm') 

Model outputs

for May only

Curvilinear

      m = lm(VCL ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'May']) 
      summary(m)
## 
## Call:
## lm(formula = VCL ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "May"])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0759  -5.7425  -0.3323   3.7169  17.8986 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        23.531      6.017   3.911 0.000348 ***
## log(motileCount)    5.309      1.143   4.643 3.67e-05 ***
## MorphIndependent    3.033      3.355   0.904 0.371452    
## MorphSatellite     -0.628      3.715  -0.169 0.866605    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.408 on 40 degrees of freedom
## Multiple R-squared:  0.3659, Adjusted R-squared:  0.3184 
## F-statistic: 7.694 on 3 and 40 DF,  p-value: 0.000357

Average path

      m = lm(VAP ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'May']) 
      summary(m)
## 
## Call:
## lm(formula = VAP ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "May"])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0345 -2.1294 -0.6012  2.8811 10.0091 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        9.9635     3.0972   3.217  0.00257 ** 
## log(motileCount)   3.8340     0.5886   6.514 8.97e-08 ***
## MorphIndependent   1.3968     1.7272   0.809  0.42345    
## MorphSatellite     0.7593     1.9121   0.397  0.69341    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.813 on 40 degrees of freedom
## Multiple R-squared:  0.5201, Adjusted R-squared:  0.4841 
## F-statistic: 14.45 on 3 and 40 DF,  p-value: 1.592e-06

Straight line

      m = lm(VSL ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'May']) 
      summary(m)
## 
## Call:
## lm(formula = VSL ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "May"])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0310 -2.4295 -0.0555  1.9077  7.0382 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        6.2035     2.8841   2.151   0.0376 *  
## log(motileCount)   3.0329     0.5481   5.534 2.14e-06 ***
## MorphIndependent   1.5611     1.6083   0.971   0.3376    
## MorphSatellite     1.8087     1.7805   1.016   0.3158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.551 on 40 degrees of freedom
## Multiple R-squared:  0.4498, Adjusted R-squared:  0.4085 
## F-statistic:  10.9 on 3 and 40 DF,  p-value: 2.295e-05
      #summary(glht(m)) 
      #plot(allEffects(m))

for June only

Curvilinear

      m = lm(VCL ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'June']) 
      summary(m)
## 
## Call:
## lm(formula = VCL ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "June"])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5518  -2.4595  -0.1141   2.5544   9.5492 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       29.4029     3.8692   7.599 3.82e-11 ***
## log(motileCount)   4.6699     0.6533   7.148 2.98e-10 ***
## MorphIndependent   0.7091     1.6673   0.425    0.672    
## MorphSatellite     0.3187     1.8025   0.177    0.860    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.398 on 84 degrees of freedom
## Multiple R-squared:  0.3799, Adjusted R-squared:  0.3578 
## F-statistic: 17.16 on 3 and 84 DF,  p-value: 8.873e-09

Average path

      m = lm(VAP ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'June']) 
      summary(m)
## 
## Call:
## lm(formula = VAP ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "June"])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4803 -1.9473  0.1188  1.7903  5.9439 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       12.2213     2.4218   5.046 2.57e-06 ***
## log(motileCount)   3.1716     0.4089   7.757 1.86e-11 ***
## MorphIndependent   1.4608     1.0436   1.400   0.1653    
## MorphSatellite     1.8965     1.1282   1.681   0.0965 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.753 on 84 degrees of freedom
## Multiple R-squared:  0.4214, Adjusted R-squared:  0.4008 
## F-statistic:  20.4 on 3 and 84 DF,  p-value: 5.081e-10

Straight line

      m = lm(VSL ~ log(motileCount) + Morph, d[Morph!='Zebra finch' & date == 'June']) 
      summary(m)
## 
## Call:
## lm(formula = VSL ~ log(motileCount) + Morph, data = d[Morph != 
##     "Zebra finch" & date == "June"])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7712 -2.5048 -0.2811  1.7716  7.8787 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        5.1769     2.8825   1.796   0.0761 .  
## log(motileCount)   2.7910     0.4867   5.735  1.5e-07 ***
## MorphIndependent   2.4964     1.2421   2.010   0.0477 *  
## MorphSatellite     3.2561     1.3428   2.425   0.0175 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.277 on 84 degrees of freedom
## Multiple R-squared:  0.3036, Adjusted R-squared:  0.2788 
## F-statistic: 12.21 on 3 and 84 DF,  p-value: 1.046e-06

for May & June

Curvilinear

      #d$Morph_F = as.factor(d$Morph)
      m = lmer(VCL ~ log(motileCount) + date + age + Morph + (1|ID), d[Morph!='Zebra finch']) 
      summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: VCL ~ log(motileCount) + date + age + Morph + (1 | ID)
##    Data: d[Morph != "Zebra finch"]
## 
## REML criterion at convergence: 810.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.07049 -0.43441 -0.04127  0.43532  2.96482 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept)  8.586   2.930   
##  Residual             22.051   4.696   
## Number of obs: 132, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       30.4545     3.6269   8.397
## log(motileCount)   4.5066     0.5845   7.710
## dateMay           -1.6869     0.9884  -1.707
## age               -0.1568     0.1926  -0.814
## MorphIndependent   1.4503     1.7713   0.819
## MorphSatellite     0.1022     1.9256   0.053
## 
## Correlation of Fixed Effects:
##             (Intr) lg(mC) dateMy age    MrphIn
## log(mtlCnt) -0.858                            
## dateMay     -0.342  0.333                     
## age         -0.228  0.009 -0.177              
## MrphIndpndn -0.455  0.038  0.074 -0.018       
## MorphSatllt -0.412  0.013  0.043  0.045  0.793

Average path

      m = lmer(VAP ~ log(motileCount) + date + age + Morph + (1|ID), d[Morph!='Zebra finch'])
      summary(m) 
## Linear mixed model fit by REML ['lmerMod']
## Formula: VAP ~ log(motileCount) + date + age + Morph + (1 | ID)
##    Data: d[Morph != "Zebra finch"]
## 
## REML criterion at convergence: 668.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8749 -0.6620  0.0003  0.5361  3.1692 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1.339    1.157   
##  Residual             8.361    2.892   
## Number of obs: 132, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       11.5733     2.0348   5.688
## log(motileCount)   3.4044     0.3326  10.235
## dateMay            0.7557     0.5933   1.274
## age               -0.1138     0.1055  -1.079
## MorphIndependent   1.4643     0.9519   1.538
## MorphSatellite     1.4853     1.0349   1.435
## 
## Correlation of Fixed Effects:
##             (Intr) lg(mC) dateMy age    MrphIn
## log(mtlCnt) -0.869                            
## dateMay     -0.347  0.326                     
## age         -0.218  0.008 -0.191              
## MrphIndpndn -0.432  0.036  0.079 -0.027       
## MorphSatllt -0.394  0.015  0.048  0.035  0.788

Straight line

      m = lmer(VSL ~ log(motileCount) + date + age + Morph + (1|ID), d[Morph!='Zebra finch']) 
      summary(m)
## Linear mixed model fit by REML ['lmerMod']
## Formula: VSL ~ log(motileCount) + date + age + Morph + (1 | ID)
##    Data: d[Morph != "Zebra finch"]
## 
## REML criterion at convergence: 686.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7980 -0.6785 -0.1253  0.5520  2.1217 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 1.811    1.346   
##  Residual             9.414    3.068   
## Number of obs: 132, groups:  ID, 91
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)       5.39488    2.19156   2.462
## log(motileCount)  2.87049    0.35750   8.029
## dateMay           1.33329    0.63226   2.109
## age              -0.05659    0.11403  -0.496
## MorphIndependent  2.12757    1.03220   2.061
## MorphSatellite    2.62987    1.12220   2.343
## 
## Correlation of Fixed Effects:
##             (Intr) lg(mC) dateMy age    MrphIn
## log(mtlCnt) -0.868                            
## dateMay     -0.347  0.327                     
## age         -0.219  0.008 -0.189              
## MrphIndpndn -0.436  0.037  0.078 -0.026       
## MorphSatllt -0.397  0.015  0.047  0.037  0.789
      #m = lm(VSL ~ log(motileCount) + date + age + Morph, d[Morph!='Zebra finch']) 
      #summary(glht(m)) 
      #summary(m)
      #plot(allEffects(m))


Test Morph-blind

    densityplot(d$VCL)
    d[VCL<40]
    sample(c('a','b','c'), size = nrow(d), replace = TRUE)
    d[, blind := sample(c('a','b','c'), size = nrow(d), replace = TRUE)]
    ggplot(d,aes(x = blind, y = VCL)) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), col = 'grey', aes(fill =blind), dotsize = 2)+
      geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
      stat_summary(fun=mean, geom="point", color="red", fill="red") +
      scale_fill_viridis(discrete=TRUE)+
      ylab('Curvilinear velocity') +
      theme_bw() +
      #guides(x =  guide_axis(angle = -45)) +
      theme(legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )

    m = lm(VCL ~ log(motileCount) + blind, d) 
    summary(glht(m)) 

# END